Exploring metabolism-growth relationships within and across individual fish.
E. Hoots and P. Biro
2025-06-09
READ ME
This R code is used to take raw respirometry slopes measured on Galaxias maculatus at Deakin’s Marine Research Centre in Queenscliff, Victoria, Australia and calculate their standard, routine, and maximum metabolic rate. These metabolic measurements are made repeatedly on known individuals held under one of two temperature treatments (18°C or 23°C), and paired with repeated mass and length measurements over five month of data collection (May - September, 2023). After reading in all data, a multivariate mixed effects modeling approach is undertaken in order to assess covariance between growth and metabolism traits among- and within-individuals.
CONTACT
Elizabeth C. Hoots
Email: beth.hoots@research.deakin.edu.au
GitHub: https://github.com/Bhoots99
LOAD REQUIRED PACKAGES
# installs and loads packages, need to install "pacman" first
pacman::p_load("brms",
"broom",
"broom.mixed",
"car",
"coda",
"dplyr",
"ggeffects",
"ggiraphExtra",
"ggplot2",
"ggpubr",
"hms",
"lme4",
"lmerTest",
"lubridate",
"nlme",
"posterior",
"readxl",
"tidyverse"
)Working directories
wd <- getwd() # getwd tells us what the current wd is, we are using this to drop it in a variable called wd
raw_data_wd <- paste0(wd, "./00_raw_data") # creates a variable with the name of the wd where raw data are stored
MRcalcs_wd <- paste0(raw_data_wd, "./MRcalcs")
MR_slopes_wd <- paste0(raw_data_wd, "./MR_Slopes")
MMR_wd <- paste0(raw_data_wd, "./MMR")
metadata_wd <- paste0(wd, "./01_metadata") # creates a variable with the name of the wd where metadata are stored
figures_wd <- paste0(wd, "./03_figures") # creates a variable with the name of the wd where figures are to be storedRead in data, calculate metabolic rates, merge into master data table
#calcSMR function adapted from Chabot, 2016 (https://doi.org/10.1111/jfb.12845)
#note that we have not used the mlnd function as the mean of the lowest 10th percentile was better suited to our data
calcSMR <- function(Y) {
u <- sort(Y)
tenpc <- round(0.1 * length(u))
SD10pc <- sd(u[1:tenpc])
low10pc = mean(u[(which((u > (mean(u[1:tenpc])-SD10pc)))):(tenpc+which((u > (mean(u[1:tenpc])-SD10pc -u[1]))))])
return(list(low10pc = low10pc))
}
# Define a function to extract date strings from file names
extract_date <- function(files) {
gsub(".*?([0-9]{8}).*", "\\1", basename(files))
}
# Specify the directories where your files are located
csv_folder <- MR_slopes_wd
xlsx_folder <- MMR_wd
# Get the list of files in each folder
csv_files <- list.files(csv_folder, pattern = "_MR_slopes.csv", full.names = TRUE)
xlsx_files <- list.files(xlsx_folder, pattern = "_MMR.xlsx", full.names = TRUE)
# Extract date strings from file names
csv_dates <- extract_date(csv_files)
xlsx_dates <- extract_date(xlsx_files)
matching_files <- Map(function(date) {
csv_file <- csv_files[grep(date, csv_dates)]
xlsx_file <- xlsx_files[grep(date, xlsx_dates)]
return(list(csv = csv_file, xlsx = xlsx_file))
}, unique(c(csv_dates, xlsx_dates)))
# Loop through matching files and perform your data processing
for (i in seq_along(matching_files)) {
date <- unique(csv_dates)[i]
current_files <- matching_files[[i]] # Renamed 'files' to 'current_files'
# Process CSV data
tb_respirometry <- read_csv(current_files$csv) %>%
rename(
chamber_ch = Chamber.No,
ID_fish = Ind,
mass = Mass,
length = Length,
volume_ch = Ch.Volume,
DOunit = DO.unit,
dateTime = Date.Time,
Temp_class = Temp.class,
slope_wBR = Slope.with.BR,
BRSlope = BR.Slope
) %>%
mutate(
dateTime = as.POSIXct(ymd_hms(dateTime)),
Time = as_hms(ymd_hms(dateTime)),
Date = as.Date(dateTime, format = "%Y/%m/%d"),
)
tb_respirometry <- tb_respirometry %>%
mutate(
volume_net = volume_ch - mass,
MR_wBR = abs(slope_wBR)*(volume_net/1000)*60*60, #uncomment for mgO2/hr instead
BR = BRSlope*(volume_ch/1000)*60*60,
MR = MR_wBR + BR,
BR_perc = (BR/MR_wBR)*100
)
#Calculate RMR and variance in MR for each fish, create a new table for individual RMR calcs
tb_rmr <- tb_respirometry %>%
group_by(chamber_ch) %>%
arrange(chamber_ch, dateTime) %>%
slice(3:(n() - 1)) %>%
ungroup() %>%
group_by(
ID_fish, mass, length, volume_net, chamber_ch, Temp_class
) %>%
arrange(ID_fish) %>%
drop_na() %>%
summarise(
RMR = mean(MR),
RMR_var = var(MR),
RMR_perc = (RMR_var/RMR)*100,
time_start = dateTime %>% min(),
time_end = dateTime %>% max()
)
#Calculate SMR for each fish, create a new table for individual SMR calcs
tb_smr <-
tb_respirometry %>%
group_by(
ID_fish, mass,length, volume_net, chamber_ch
) %>%
arrange(ID_fish) %>%
drop_na() %>%
summarise(
SMR = calcSMR(MR)$low10pc %>% unname(),
time_start = dateTime %>% min(),
time_end = dateTime %>% max()
)
# Process XLSX data
tb_mmr <- read_excel(current_files$xlsx) %>%
rename(
chamber_ch = Chamber.No,
ID_fish = Ind,
mass = Mass,
volume_ch = Ch.Volume,
DOunit = DO.unit,
dateTime = Date.Time,
Temp_class = Temp.class,
slope_wBR = Slope.with.BR,
BRSlope = BR.Slope
) %>%
mutate(
dateTime = as.POSIXct(ymd_hms(dateTime)),
Time = as_hms(ymd_hms(dateTime)),
Date = as.Date(dateTime, format = "%Y/%m/%d"),
) %>%
drop_na()
tb_mmr <- tb_mmr %>%
mutate(
volume_net = volume_ch - mass,
MMR_wBR = abs(slope_wBR)*(volume_net/1000)*60*60, #uncomment for mgO2/hr instead
BR = abs(BRSlope)*(volume_ch/1000)*60*60,
MMR = MMR_wBR - BR,
BR_perc = (BR/MMR_wBR)*100
) %>%
arrange(ID_fish)
# Combine the processed data
# Combine the processed data, selecting only specific columns
tb_MR_master <- tb_rmr %>%
left_join(select(tb_smr, ID_fish, SMR), by = "ID_fish") %>%
left_join(select(tb_mmr, ID_fish, MMR), by = "ID_fish") %>%
select(-matches(".x$")) %>%
rename_with(~gsub("\\.y$", "", .), matches(".y$"))
# Export data
setwd(MRcalcs_wd)
write.csv(tb_MR_master, file = paste0(date, "_MRcalcs.csv"), col.names = NA, row.names = FALSE)
}
# read all "_MRcalcs.csv" files, bind into one table, and add columns with logged values
file.list <- list.files(MRcalcs_wd)
setwd(MRcalcs_wd)
calcs_all <- lapply(file.list, read_csv)
tb_MRcalcs <- bind_rows(calcs_all, .id="Resp_Day") %>%
mutate(
log_SMR_low10pc = log10(SMR),
log_mass = log10(mass),
log_RMR = log10(RMR),
log_MMR = log10(MMR),
Resp_Day = as.numeric(Resp_Day)
)
# assign Month category to all timepoints
tb_MRcalcs$Month <- cut(
tb_MRcalcs$Resp_Day,
breaks = c(0, 13, seq(23, max(tb_MRcalcs$Resp_Day) + 10, by = 10)),
right = FALSE,
labels = c("April", "May", "June", "July", "August", "September")
)
setwd(raw_data_wd)
biometrics <- read.csv("FinalBiometrics.csv")
ds <- left_join(tb_MRcalcs, biometrics, by = c("ID_fish")) #a table with all resp data and metadata collected in dissection
ds1 <- ds[which(ds$Month != "April"),] # remove April data from analysis due to temperature inconsistencyMultivariate mixed-effects models to test the relationships between growth and metabolism in Galaxias maculatus.
ds1 <- ds1 %>%
# Filter to keep only ID_fish where a mass exists for Month == "May"
filter(ID_fish %in% ID_fish[Month == "May" & !is.na(mass)])
month_order <- c("May", "June", "July", "August", "September") # assign chronological order to Month (otherwise alphabetical)
ds1$Month <- factor(ds1$Month, levels = month_order) # assign month_order to factor variable "Month" in ds1
min_mass = min(ds1$mass)
min_log_mass = min(ds1$log_mass)
min_finalmass = min(((ds1[which(!is.na(ds1$Final_Mass_g)),]$Final_Mass_g)))
min_log_finalmass = min((log10(ds1[which(!is.na(ds1$Final_Mass_g)),]$Final_Mass_g)))
# format data columns and generate centred and log-scaled variables
ds1 <- ds1 %>%
mutate(
mass_centre = mass - min_mass, # centred mass
log_mass_centre = log_mass - log10(min_mass), # log-scaled centred mass
Temp_class = as.factor(Temp_class), # set temperature treatments as factors
SMR1 = log10(SMR), # log-scaled SMR
RMR1 = log10(RMR), # log-scaled RMR
MMR1 = log10(MMR), # log-scaled MMR
mass1 = log_mass, # log-scaled mass
Month_int = case_when(Month == "May" ~ 1,
Month == "June" ~ 2,
Month == "July" ~ 3,
Month == "August" ~ 4,
Month == "September" ~ 5),
Month_centred = Month_int - 1, # left-centred mass so May is the intercept
Month_centred_factor = as.factor(Month_centred), # set centred month as factors
Rack = as.integer(Rack), # set rack ID to integer (for now)
Bin = as.integer(Bin), # set bin ID to integer (for now)
finalmass_centre = Final_Mass_g - min_finalmass, # centre last mass measurement taken
log_finalmass_centre = log10(Final_Mass_g) - min_log_finalmass, # log-scaled and centred final mass
log_fatmass = log10(Fat_mass_g) # log-scaled fat mass
)
ds1 <- tidyr::unite(data = ds1, Rack, Bin, col = "tankID", sep = "_") # merge rack and bin data into one identifier, call it tankID
#read in density data table, call it BinCounts
setwd(raw_data_wd)
BinCounts <- read_excel("BinCounts.xlsx", sheet = "data")
#merge BinCounts and ds1
ds1 <- left_join(ds1, BinCounts %>% dplyr::select("Month", "tankID", "Population"), by = c("Month", "tankID"))
#reformat Sex variable so that both unknown and immature fish are "NA"
ds1$Sex <- ifelse(ds1$Sex == "I", NA, ds1$Sex)
# rename "population" as density, set tankID as factor
ds1 <- ds1 %>%
rename(
Density = Population
) %>%
mutate(
tankID = as.factor(tankID)
)
# check structure of ds1
str(ds1)## tibble [624 × 41] (S3: tbl_df/tbl/data.frame)
## $ Resp_Day : num [1:624] 13 13 13 13 13 13 13 13 13 13 ...
## $ ID_fish : num [1:624] 81 82 83 84 85 86 87 88 89 90 ...
## $ chamber_ch : chr [1:624] "A1" "B3" "B4" "A2" ...
## $ Temp_class : Factor w/ 2 levels "18","23": 2 2 2 2 2 2 2 2 2 2 ...
## $ RMR : num [1:624] 0.952 0.412 0.712 1.016 0.876 ...
## $ RMR_var : num [1:624] 0.0149 0.0135 0.0429 0.0696 0.0454 ...
## $ RMR_perc : num [1:624] 1.57 3.28 6.02 6.85 5.18 ...
## $ time_start : POSIXct[1:624], format: "2023-05-22 17:53:00" "2023-05-22 17:53:00" ...
## $ time_end : POSIXct[1:624], format: "2023-05-23 08:03:00" "2023-05-23 08:03:00" ...
## $ mass : num [1:624] 2.65 0.81 1.78 2.73 2.29 0.98 2.04 2.43 2.31 1.65 ...
## $ length : num [1:624] 77 62 73 78 80 64 80 77 78 68 ...
## $ volume_net : num [1:624] 707 709 708 707 708 ...
## $ SMR : num [1:624] 0.845 0.284 0.53 0.804 0.707 ...
## $ MMR : num [1:624] 1.43 1.11 1.44 1.64 1.69 ...
## $ log_SMR_low10pc : num [1:624] -0.0733 -0.5461 -0.2758 -0.0948 -0.1506 ...
## $ log_mass : num [1:624] 0.4232 -0.0915 0.2504 0.4362 0.3598 ...
## $ log_RMR : num [1:624] -0.02124 -0.38552 -0.14722 0.00674 -0.05726 ...
## $ log_MMR : num [1:624] 0.1544 0.0458 0.1591 0.2154 0.2275 ...
## $ Month : chr [1:624] "May" "May" "May" "May" ...
## $ tankID : Factor w/ 20 levels "1_10","1_11",..: 17 17 17 17 17 17 17 17 18 18 ...
## $ TagID : chr [1:624] "YL" "YR" "BrL" "BrR" ...
## $ Final_Mass_g : num [1:624] 4.62 3.65 2.64 6.33 5.19 ...
## $ Final_Length_mm : int [1:624] 94 86 84 98 100 95 90 89 91 98 ...
## $ Gonad_mass_g : num [1:624] 0.035 0.022 0.018 0.039 0.052 0.028 0.362 0.033 0.015 0.093 ...
## $ Sex : chr [1:624] "F" "M" "F" "M" ...
## $ Fat_mass_g : num [1:624] 0.315 0.123 0.062 0.314 0.214 0.091 0.076 0.044 0.15 0.209 ...
## $ Ventricle_mass_g : chr [1:624] "0.017" "" "" "" ...
## $ Otolith_age_y : num [1:624] 1 NA NA NA NA 1 1 NA NA 1 ...
## $ mass_centre : num [1:624] 2.07 0.23 1.2 2.15 1.71 0.4 1.46 1.85 1.73 1.07 ...
## $ log_mass_centre : num [1:624] 0.66 0.145 0.487 0.673 0.596 ...
## $ SMR1 : num [1:624] -0.0733 -0.5461 -0.2758 -0.0948 -0.1506 ...
## $ RMR1 : num [1:624] -0.02124 -0.38552 -0.14722 0.00674 -0.05726 ...
## $ MMR1 : num [1:624] 0.1544 0.0458 0.1591 0.2154 0.2275 ...
## $ mass1 : num [1:624] 0.4232 -0.0915 0.2504 0.4362 0.3598 ...
## $ Month_int : num [1:624] 1 1 1 1 1 1 1 1 1 1 ...
## $ Month_centred : num [1:624] 0 0 0 0 0 0 0 0 0 0 ...
## $ Month_centred_factor: Factor w/ 5 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ finalmass_centre : num [1:624] 3.3 2.33 1.32 5.01 3.87 ...
## $ log_finalmass_centre: num [1:624] 0.544 0.442 0.301 0.681 0.594 ...
## $ log_fatmass : num [1:624] -0.502 -0.91 -1.208 -0.503 -0.67 ...
## $ Density : num [1:624] 8 8 8 8 8 8 8 8 5 5 ...head(ds1)## # A tibble: 6 × 41
## Resp_Day ID_fish chamber_ch Temp_class RMR RMR_var RMR_perc
## <dbl> <dbl> <chr> <fct> <dbl> <dbl> <dbl>
## 1 13 81 A1 23 0.952 0.0149 1.57
## 2 13 82 B3 23 0.412 0.0135 3.28
## 3 13 83 B4 23 0.712 0.0429 6.02
## 4 13 84 A2 23 1.02 0.0696 6.85
## 5 13 85 A4 23 0.876 0.0454 5.18
## 6 13 86 A3 23 0.751 0.0504 6.71
## # ℹ 34 more variables: time_start <dttm>, time_end <dttm>, mass <dbl>,
## # length <dbl>, volume_net <dbl>, SMR <dbl>, MMR <dbl>,
## # log_SMR_low10pc <dbl>, log_mass <dbl>, log_RMR <dbl>, log_MMR <dbl>,
## # Month <chr>, tankID <fct>, TagID <chr>, Final_Mass_g <dbl>,
## # Final_Length_mm <int>, Gonad_mass_g <dbl>, Sex <chr>, Fat_mass_g <dbl>,
## # Ventricle_mass_g <chr>, Otolith_age_y <dbl>, mass_centre <dbl>,
## # log_mass_centre <dbl>, SMR1 <dbl>, RMR1 <dbl>, MMR1 <dbl>, mass1 <dbl>, …GROWTH ANALYSES
Visualize growth data
hist(ds1$mass) # plot all mass values to see mass range and most frequent values
ds1$Month <- factor(ds1$Month, levels = month_order) # reassign month_order to factor variable "Month" in ds1
# Plot mass over time, connect by ID_fish to show individual trends
ggplot(ds1, aes(x = Month, y = mass, group = ID_fish)) + # plot raw data by ID
geom_point(size = 2) + # size of data points
geom_line() + # connect points by individual ID
facet_wrap(. ~ Temp_class) + # produces one plot per temperature treatment
theme_bw()
Log-linearize mass data
We can also visualise this by plotting log10 mass over log10_x (note that we cannot use log10_y as some masses were < 1 initially)
ggplot(ds1, aes(x = Month_int, y = log_mass, group = ID_fish)) + # plot logged mass data by ID
geom_point(size = 2) + # size of data points
geom_line() + # connect points by individual ID
facet_wrap(. ~ Temp_class) + # produces one plot per temperature treatment
theme_bw() +
scale_x_log10() # display x-axis on log scale
Univariate growth models for two temperature treatments
ds18 <- ds1[which(ds1$Temp_class == "18"),] # subset ds for only fish at 18C
ds23 <- ds1[which(ds1$Temp_class == "23"),] # subset ds for only fish at 23C
mass18_mod<-lme4::lmer(log_mass ~ Month + (1 + Month_int|ID_fish), # model structure, 18C only
data=ds18, REML = T)
summary(mass18_mod) # generate model output## Linear mixed model fit by REML ['lmerMod']
## Formula: log_mass ~ Month + (1 + Month_int | ID_fish)
## Data: ds18
##
## REML criterion at convergence: -492
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8214 -0.4648 0.0139 0.4865 2.3450
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 0.040770 0.20191
## Month_int 0.002038 0.04515 -0.57
## Residual 0.003528 0.05939
## Number of obs: 317, groups: ID_fish, 72
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.32789 0.02232 14.689
## MonthJune 0.15461 0.01147 13.478
## MonthJuly 0.22113 0.01498 14.761
## MonthAugust 0.15515 0.01989 7.799
## MonthSeptember 0.36370 0.02487 14.622
##
## Correlation of Fixed Effects:
## (Intr) MnthJn MnthJl MnthAg
## MonthJune -0.364
## MonthJuly -0.412 0.629
## MonthAugust -0.410 0.600 0.770
## MonthSptmbr -0.408 0.581 0.776 0.865plot(mass18_mod) # residuals plot
mass23_mod<-lme4::lmer(log_mass ~ Month + (1 + Month_int|ID_fish), # model structure, 23C only
data=ds23, REML=T)
summary(mass23_mod) # generate model output## Linear mixed model fit by REML ['lmerMod']
## Formula: log_mass ~ Month + (1 + Month_int | ID_fish)
## Data: ds23
##
## REML criterion at convergence: -515.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9289 -0.3941 -0.0018 0.5057 4.8187
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 0.040995 0.20247
## Month_int 0.001445 0.03802 -0.46
## Residual 0.003138 0.05602
## Number of obs: 307, groups: ID_fish, 64
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.21996 0.02451 8.974
## MonthJune 0.12091 0.01123 10.771
## MonthJuly 0.21448 0.01402 15.297
## MonthAugust 0.22778 0.01773 12.848
## MonthSeptember 0.35157 0.02193 16.028
##
## Correlation of Fixed Effects:
## (Intr) MnthJn MnthJl MnthAg
## MonthJune -0.298
## MonthJuly -0.335 0.623
## MonthAugust -0.341 0.611 0.773
## MonthSptmbr -0.337 0.589 0.778 0.858plot(mass23_mod) # residuals plot
Combined growth model
# mass model including both temperature treatments, and their interaction with Month
mass_mod <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),
data=ds1, REML=T)
summary(mass_mod) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (1 +
## Month_int | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -789
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6831 -0.4303 0.0066 0.4509 4.5846
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 0.038839 0.19708
## Month_int 0.001745 0.04177 -0.55
## Residual 0.003485 0.05903
## Number of obs: 509, groups: ID_fish, 113
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.175461 0.059281 337.150651 2.960
## MonthJune 0.156983 0.011479 391.676251 13.676
## MonthJuly 0.224031 0.014800 247.773510 15.137
## MonthAugust 0.162564 0.020135 174.402058 8.074
## MonthSeptember 0.373693 0.024476 123.702938 15.268
## as.factor(Temp_class)23 -0.071983 0.036905 120.004217 -1.950
## SexF 0.123289 0.043645 115.959487 2.825
## SexM 0.051217 0.046603 115.166678 1.099
## Density 0.010421 0.005792 468.638813 1.799
## MonthJune:as.factor(Temp_class)23 -0.026407 0.018735 390.916591 -1.409
## MonthJuly:as.factor(Temp_class)23 -0.003931 0.023841 236.415811 -0.165
## MonthAugust:as.factor(Temp_class)23 0.065838 0.031139 146.333633 2.114
## MonthSeptember:as.factor(Temp_class)23 -0.015746 0.038370 107.376799 -0.410
## Pr(>|t|)
## (Intercept) 0.00330 **
## MonthJune < 2e-16 ***
## MonthJuly < 2e-16 ***
## MonthAugust 1.07e-13 ***
## MonthSeptember < 2e-16 ***
## as.factor(Temp_class)23 0.05345 .
## SexF 0.00557 **
## SexM 0.27405
## Density 0.07265 .
## MonthJune:as.factor(Temp_class)23 0.15948
## MonthJuly:as.factor(Temp_class)23 0.86916
## MonthAugust:as.factor(Temp_class)23 0.03618 *
## MonthSeptember:as.factor(Temp_class)23 0.68236
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(mass_mod) # residuals plot
ranova(mass_mod) # ANOVA table output to test random effects## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + Month:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 17 394.50 -755.01
## Month_int in (1 + Month_int | ID_fish) 15 331.56 -633.11 125.9 2 < 2.2e-16
##
## <none>
## Month_int in (1 + Month_int | ID_fish) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# mass model including both temperature treatments, and their interaction with Month
mass_mod_noMonth <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex + Density + (1|ID_fish),
data=ds1, REML=T)
summary(mass_mod_noMonth) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (1 |
## ID_fish)
## Data: ds1
##
## REML criterion at convergence: -663.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6703 -0.5089 0.0290 0.4797 3.9785
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.026454 0.16265
## Residual 0.007503 0.08662
## Number of obs: 509, groups: ID_fish, 113
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.193236 0.060932 357.632744 3.171
## MonthJune 0.156704 0.015013 391.702422 10.438
## MonthJuly 0.221347 0.015486 395.495181 14.293
## MonthAugust 0.161419 0.017561 412.738321 9.192
## MonthSeptember 0.373567 0.017561 412.738321 21.273
## as.factor(Temp_class)23 -0.073410 0.036398 155.881117 -2.017
## SexF 0.121574 0.043722 118.677173 2.781
## SexM 0.050315 0.046686 116.817268 1.078
## Density 0.008195 0.006125 491.655639 1.338
## MonthJune:as.factor(Temp_class)23 -0.026128 0.024591 391.220473 -1.062
## MonthJuly:as.factor(Temp_class)23 -0.001748 0.024840 392.469865 -0.070
## MonthAugust:as.factor(Temp_class)23 0.065981 0.026055 399.871498 2.532
## MonthSeptember:as.factor(Temp_class)23 -0.015885 0.026056 398.802053 -0.610
## Pr(>|t|)
## (Intercept) 0.00165 **
## MonthJune < 2e-16 ***
## MonthJuly < 2e-16 ***
## MonthAugust < 2e-16 ***
## MonthSeptember < 2e-16 ***
## as.factor(Temp_class)23 0.04543 *
## SexF 0.00631 **
## SexM 0.28337
## Density 0.18154
## MonthJune:as.factor(Temp_class)23 0.28868
## MonthJuly:as.factor(Temp_class)23 0.94394
## MonthAugust:as.factor(Temp_class)23 0.01171 *
## MonthSeptember:as.factor(Temp_class)23 0.54244
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(mass_mod_noMonth) # residuals plot
ranova(mass_mod_noMonth) # ANOVA table output to test random effects## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + Month:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 15 331.56 -633.11
## (1 | ID_fish) 14 115.72 -203.43 431.68 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# mass model with random intercepts by temp treatment
mass_mod_trtint <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex + Density + (0 + Month_int|Temp_class),
data=ds1, REML=T)
summary(mass_mod_trtint) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (0 +
## Month_int | Temp_class)
## Data: ds1
##
## REML criterion at convergence: -231.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.77281 -0.70174 0.00548 0.69951 2.82668
##
## Random effects:
## Groups Name Variance Std.Dev.
## Temp_class Month_int 0.12294 0.3506
## Residual 0.03292 0.1814
## Number of obs: 509, groups: Temp_class, 2
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 2.277e-01 3.554e-01 4.960e+02 0.641
## MonthJune 1.658e-01 3.520e-01 4.960e+02 0.471
## MonthJuly 2.258e-01 7.020e-01 4.960e+02 0.322
## MonthAugust 1.583e-01 1.052e+00 4.960e+02 0.150
## MonthSeptember 3.704e-01 1.403e+00 4.960e+02 0.264
## as.factor(Temp_class)23 -7.388e-02 4.971e-01 4.960e+02 -0.149
## SexF 1.055e-01 2.758e-02 4.960e+02 3.824
## SexM 3.331e-02 2.858e-02 4.960e+02 1.165
## Density 5.294e-03 6.385e-03 4.960e+02 0.829
## MonthJune:as.factor(Temp_class)23 -2.945e-02 4.985e-01 4.960e+02 -0.059
## MonthJuly:as.factor(Temp_class)23 -4.927e-04 9.930e-01 4.960e+02 0.000
## MonthAugust:as.factor(Temp_class)23 7.486e-02 1.489e+00 4.960e+02 0.050
## MonthSeptember:as.factor(Temp_class)23 -1.122e-02 1.984e+00 4.960e+02 -0.006
## Pr(>|t|)
## (Intercept) 0.521955
## MonthJune 0.637892
## MonthJuly 0.747808
## MonthAugust 0.880518
## MonthSeptember 0.791861
## as.factor(Temp_class)23 0.881916
## SexF 0.000148 ***
## SexM 0.244381
## Density 0.407426
## MonthJune:as.factor(Temp_class)23 0.952906
## MonthJuly:as.factor(Temp_class)23 0.999604
## MonthAugust:as.factor(Temp_class)23 0.959909
## MonthSeptember:as.factor(Temp_class)23 0.995491
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(mass_mod_trtint) # residuals plot
ranova(mass_mod_trtint) # ANOVA table output to test random effects## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (0 + Month_int | Temp_class) + Month:as.factor(Temp_class)
## npar logLik AIC LRT Df
## <none> 15 115.72 -201.43
## Month_int in (0 + Month_int | Temp_class) 15 115.72 -201.43 5.0591e-12 0
## Pr(>Chisq)
## <none>
## Month_int in (0 + Month_int | Temp_class)Plot predicted mean level trend in mass over time, by temperature
pred <- ggpredict(mass_mod, terms = c("Month", "Temp_class")) # generate trend lines for the two temperature treatments
ggplot() +
geom_errorbar(
data = pred,
aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) + # add 95% confidence interval bars
geom_point(
data = pred,
aes(x = x, y = predicted, color = group, group = group)) + # predicted level by temp*month
geom_line(
data = pred,
aes(x = x, y = predicted, color = group, group = group)) + # show level trend over time
labs(
x = "Month",
y = "Predicted log_mass",
color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00"))
F test on interaction effect
anova(mass_mod, type = "3") # type 3 anova because we expect a meaningful interaction## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Month 1.80879 0.45220 4 215.74 129.7554 < 2.2e-16 ***
## as.factor(Temp_class) 0.01454 0.01454 1 107.54 4.1733 0.04351 *
## Sex 0.03252 0.01626 2 110.81 4.6664 0.01133 *
## Density 0.01128 0.01128 1 468.64 3.2366 0.07265 .
## Month:as.factor(Temp_class) 0.10964 0.02741 4 212.06 7.8649 6.274e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Generate fitted data
z_mass<-augment(mass_mod) # get variables and fitted values from model
head(z_mass)## # A tibble: 6 × 19
## .rownames log_mass Month `as.factor(Temp_class)` Sex Density Month_int
## <chr> <dbl> <fct> <fct> <fct> <dbl> <dbl>
## 1 1 0.423 May 23 F 8 1
## 2 2 -0.0915 May 23 M 8 1
## 3 3 0.250 May 23 F 8 1
## 4 4 0.436 May 23 M 8 1
## 5 5 0.360 May 23 M 8 1
## 6 7 0.310 May 23 M 8 1
## # ℹ 12 more variables: ID_fish <dbl>, .fitted <dbl>, .resid <dbl>, .hat <dbl>,
## # .cooksd <dbl>, .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>,
## # .sqrtrwt <dbl>, .weights <dbl>, .wtres <dbl>##### check assumptions: plot all residual values against the fitted
ggplot(z_mass, aes(x = .fitted, y = .resid)) + # same plot as before
geom_point(size = 2) 
FIGURE 1
fig1.1 <- ggplot(z_mass, aes(x = Month, y = log_mass, group = ID_fish)) +
geom_line(aes(y = .fitted, color = `as.factor(Temp_class)`), # note this will look odd if there are additional fixed effects
alpha = 0.6) +
theme_classic() +
labs(x = "Month", y = expression(log[10]~(mass~(g))), color = expression("Temperature ("*degree*C*")")) +
ylim(-0.24,1) +
scale_color_manual(values = c("#78B7C5", "#F21A00"))
fig1.2 <- ggplot() +
geom_errorbar(
data = pred,
aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2, alpha = 0.3) + # add 95% confidence interval bars
geom_line(
data = pred,
aes(x = x, y = predicted, color = group, group = group),
size = 1) +
labs(
x = "Month",
y = expression(log[10]~(mass~(g))),
color = expression("Temperature ("*degree*C*")")) +
theme_classic() +
ylim(-0.24,1) +
scale_color_manual(values = c("#78B7C5", "#F21A00"))
fig1 <- ggarrange(print(fig1.1 + rremove("xlab") + rremove ("ylab")),
print(fig1.2 + rremove("xlab") + rremove ("ylab")),
ncol = 2, nrow = 1, widths = c(1, 1.1),
align = "v",
common.legend = TRUE,
legend = "right")
require(grid)
jpeg(file = "Fig1.jpeg", width=960, height=480, units = "px")
annotate_figure(
fig1,
bottom = textGrob("Month"),
left = textGrob(
expression(log[10]~(mass~(g))),
rot = 90, # Rotate the y-axis label 90 degrees
gp = gpar(fontsize = 12) # Optionally adjust the font size
)
)
setwd(figures_wd)
jpeg(file = "Fig1.jpeg", width=960, height=480, units = "px")
annotate_figure(
fig1,
bottom = textGrob("Month"),
left = textGrob(
expression(log[10]~(mass~(g))),
rot = 90, # Rotate the y-axis label 90 degrees
gp = gpar(fontsize = 12) # Optionally adjust the font size
)
)
dev.off()## png
## 2SMR MODEL
Visualize SMR data
hist(ds1$SMR) # plot all raw SMR values in a histogram, note skewness
hist(log10(ds1$SMR)) # plot all log-SMR values in a histogram
ggplot(ds1, aes(x = mass, y = SMR, group = ID_fish)) + # plot raw SMR data by ID against mass
geom_point(size = 2) + # size of data points
facet_wrap(. ~ Temp_class) + # produces one plot per temperature treatment
theme_bw() 
Visualize log-mass and SMR relationship
ggplot(ds1, aes(x = log_mass, y = SMR, group = ID_fish)) + # plot raw SMR data by ID against log10(mass)
geom_point(size = 2) + # size of data points
facet_wrap(. ~ Temp_class) + # produces one plot per temperature treatment
theme_bw() 
Note that this relationship is non-linear, we expect this to be described by a power function y = a*x^b
This reflects our understanding of the relationship between mass and metabolism, which can be described by the power function y = a*x^b
ggplot(ds1, aes(x = log_mass, y = SMR, group = ID_fish)) + # plot raw SMR data by ID against log10(mass)
geom_point(size = 2) + # size of data points
facet_wrap(. ~ Temp_class) + # produces one plot per temperature treatment
theme_bw() +
scale_y_log10() # plot data with log10 scaled y axis
The trends become linear when SMR is presented on log scale
Univariate SMR models for two temperature treatments
# model structure using log10 SMR and mass for fish at 18C only
SMR18_mod<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish), # random slopes by month, random intercepts by ID_fish
data=ds18, REML = T)
summary(SMR18_mod) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
## Data: ds18
##
## REML criterion at convergence: -553.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6860 -0.7401 0.0462 0.7148 2.7433
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 1.063e-17 3.261e-09
## Residual 9.832e-03 9.916e-02
## Number of obs: 317, groups: ID_fish, 72
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.67878 0.01415 315.00000 -47.96 <2e-16 ***
## log_mass 0.78058 0.02562 315.00000 30.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## log_mass -0.919
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')plot(SMR18_mod) # plot model residuals
ranova(SMR18_mod) # anova-like test shows insignificant random slopes## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 4 276.78 -545.56
## (1 | ID_fish) 3 276.78 -547.56 0 1 1SMR18_mod_1<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish), # random slopes term removed
data=ds18, REML = T)
summary(SMR18_mod_1) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
## Data: ds18
##
## REML criterion at convergence: -553.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6860 -0.7401 0.0462 0.7148 2.7433
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 1.063e-17 3.261e-09
## Residual 9.832e-03 9.916e-02
## Number of obs: 317, groups: ID_fish, 72
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.67878 0.01415 315.00000 -47.96 <2e-16 ***
## log_mass 0.78058 0.02562 315.00000 30.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## log_mass -0.919
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')plot(SMR18_mod_1) # plot model residuals
ranova(SMR18_mod_1) # anova-like test shows insignificant random intercepts, but retain in model for testing## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 4 276.78 -545.56
## (1 | ID_fish) 3 276.78 -547.56 0 1 1# model structure using log10 SMR and mass for fish at 23C only
SMR23_mod<-lmer(log10(SMR) ~ log_mass + (1 + Month_int|ID_fish), # random slopes by month, random intercepts by ID_fish
data=ds23, REML=T)
summary(SMR23_mod) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 + Month_int | ID_fish)
## Data: ds23
##
## REML criterion at convergence: -649.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1396 -0.5503 -0.0023 0.5686 3.2950
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 0.0052305 0.07232
## Month_int 0.0003171 0.01781 -0.94
## Residual 0.0056041 0.07486
## Number of obs: 307, groups: ID_fish, 64
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.47671 0.01133 101.03860 -42.07 <2e-16 ***
## log_mass 0.61389 0.02319 128.47109 26.48 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## log_mass -0.876plot(SMR23_mod) # plot model residuals
ranova(SMR23_mod) # anova-like test shows insignificant random slopes## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + (1 + Month_int | ID_fish)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 6 324.70 -637.41
## Month_int in (1 + Month_int | ID_fish) 4 320.96 -633.92 7.4863 2 0.02368
##
## <none>
## Month_int in (1 + Month_int | ID_fish) *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1SMR23_mod_1<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish), # random slopes term removed
data=ds23, REML=T)
summary(SMR23_mod_1) # generate model output## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
## Data: ds23
##
## REML criterion at convergence: -641.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0810 -0.5308 -0.0445 0.5962 3.1721
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.0006479 0.02545
## Residual 0.0064034 0.08002
## Number of obs: 307, groups: ID_fish, 64
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.47203 0.01064 133.59553 -44.38 <2e-16 ***
## log_mass 0.60664 0.02244 178.85157 27.03 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## log_mass -0.851plot(SMR23_mod_1) # plot model residuals
ranova(SMR23_mod_1) # anova-like test shows insignificant random intercepts, but retain in model for testing## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 4 320.96 -633.92
## (1 | ID_fish) 3 318.78 -631.56 4.3593 1 0.03681 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Combined SMR model
# model using log10 SMR and the interaction of log10 mass and temperature treatment
SMR_mod <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|ID_fish),
data=ds1, REML=T)
plot(SMR_mod) # generate model output
summary(SMR_mod) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -927.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.82335 -0.64895 0.01139 0.65120 3.06128
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.0002204 0.01485
## Residual 0.0084623 0.09199
## Number of obs: 509, groups: ID_fish, 113
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.648757 0.030786 191.298737 -21.073
## log_mass 0.791230 0.025593 308.700979 30.916
## as.factor(Temp_class)23 0.178487 0.020589 235.734567 8.669
## SexF -0.025831 0.014616 181.956942 -1.767
## SexM -0.013551 0.014938 167.180010 -0.907
## Density -0.002235 0.003289 173.189264 -0.680
## log_mass:as.factor(Temp_class)23 -0.126302 0.039605 277.914923 -3.189
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 7.13e-16 ***
## SexF 0.07886 .
## SexM 0.36562
## Density 0.49769
## log_mass:as.factor(Temp_class)23 0.00159 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.449
## as.fc(T_)23 -0.379 0.622
## SexF -0.392 -0.194 -0.067
## SexM -0.413 -0.066 -0.014 0.790
## Density -0.834 0.120 0.120 0.123 0.109
## lg_:.(T_)23 0.340 -0.637 -0.898 0.046 -0.038 -0.105ranova(SMR_mod) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 9 463.94 -909.88
## (1 | ID_fish) 8 463.67 -911.34 0.54165 1 0.4618# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
SMR_mod_month <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),
data=ds1, REML=T)
plot(SMR_mod_month) # generate model output
summary(SMR_mod_month) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 + Month_int | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -929.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.85707 -0.62933 0.03354 0.64009 3.03111
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 1.442e-03 0.037971
## Month_int 4.593e-05 0.006777 -1.00
## Residual 8.259e-03 0.090880
## Number of obs: 509, groups: ID_fish, 113
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.644790 0.031143 232.514432 -20.704
## log_mass 0.793677 0.025854 340.474164 30.699
## as.factor(Temp_class)23 0.178273 0.021153 199.523001 8.428
## SexF -0.026653 0.015115 151.602600 -1.763
## SexM -0.015453 0.015411 150.798427 -1.003
## Density -0.002854 0.003294 238.799373 -0.866
## log_mass:as.factor(Temp_class)23 -0.128338 0.039994 319.841634 -3.209
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 6.95e-15 ***
## SexF 0.07985 .
## SexM 0.31760
## Density 0.38713
## log_mass:as.factor(Temp_class)23 0.00147 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.452
## as.fc(T_)23 -0.379 0.622
## SexF -0.404 -0.197 -0.065
## SexM -0.427 -0.069 -0.008 0.799
## Density -0.826 0.120 0.115 0.128 0.114
## lg_:.(T_)23 0.345 -0.638 -0.902 0.047 -0.040 -0.107
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')ranova(SMR_mod_month) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 464.62 -907.24
## Month_int in (1 + Month_int | ID_fish) 9 463.94 -909.88 1.3653 2 0.5053# model with random slopes by temp treatment
SMR_mod_trt <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0+Temp_class|ID_fish),
data=ds1, REML=T)
plot(SMR_mod_trt) # generate model output
summary(SMR_mod_trt) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (0 + Temp_class | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -927.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.83458 -0.62750 0.01992 0.65194 3.04447
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish Temp_class18 0.0001802 0.01342
## Temp_class23 0.0002780 0.01667 0.53
## Residual 0.0084649 0.09201
## Number of obs: 509, groups: ID_fish, 113
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.648557 0.030753 188.423517 -21.089
## log_mass 0.791050 0.025482 219.289631 31.044
## as.factor(Temp_class)23 0.178657 0.020645 208.685088 8.654
## SexF -0.025952 0.014606 182.703141 -1.777
## SexM -0.013653 0.014928 167.741231 -0.915
## Density -0.002239 0.003292 171.722340 -0.680
## log_mass:as.factor(Temp_class)23 -0.126667 0.039728 233.094125 -3.188
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 1.35e-15 ***
## SexF 0.07727 .
## SexM 0.36170
## Density 0.49730
## log_mass:as.factor(Temp_class)23 0.00163 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.448
## as.fc(T_)23 -0.374 0.617
## SexF -0.391 -0.196 -0.067
## SexM -0.413 -0.065 -0.013 0.789
## Density -0.835 0.120 0.119 0.122 0.108
## lg_:.(T_)23 0.336 -0.632 -0.898 0.047 -0.037 -0.104ranova(SMR_mod_trt) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df
## <none> 11 463.95 -905.91
## Temp_class in (0 + Temp_class | ID_fish) 9 463.94 -909.88 0.028652 2
## Pr(>Chisq)
## <none>
## Temp_class in (0 + Temp_class | ID_fish) 0.9858# model with random intercepts by temp treatment
SMR_mod_trtint <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),
data=ds1, REML=T)
plot(SMR_mod_trtint) # generate model output
summary(SMR_mod_trtint) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 | Temp_class)
## Data: ds1
##
## REML criterion at convergence: -927.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.86266 -0.65384 0.00029 0.65027 3.09859
##
## Random effects:
## Groups Name Variance Std.Dev.
## Temp_class (Intercept) 0.0002948 0.01717
## Residual 0.0086736 0.09313
## Number of obs: 509, groups: Temp_class, 2
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -6.473e-01 3.448e-02 1.288e-08 -18.773
## log_mass 7.897e-01 2.513e-02 5.020e+02 31.428
## as.factor(Temp_class)23 1.768e-01 3.152e-02 2.247e-09 5.609
## SexF -2.541e-02 1.419e-02 5.020e+02 -1.791
## SexM -1.322e-02 1.447e-02 5.020e+02 -0.913
## Density -2.402e-03 3.186e-03 5.020e+02 -0.754
## log_mass:as.factor(Temp_class)23 -1.225e-01 3.879e-02 5.020e+02 -3.158
## Pr(>|t|)
## (Intercept) 1.00000
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 1.00000
## SexF 0.07398 .
## SexM 0.36158
## Density 0.45117
## log_mass:as.factor(Temp_class)23 0.00168 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.390
## as.fc(T_)23 -0.483 0.399
## SexF -0.343 -0.194 -0.041
## SexM -0.363 -0.065 -0.006 0.797
## Density -0.720 0.114 0.076 0.126 0.111
## lg_:.(T_)23 0.298 -0.639 -0.576 0.045 -0.039 -0.105ranova(SMR_mod_trtint) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 9 463.67 -909.34
## (1 | Temp_class) 8 463.67 -911.34 1.1369e-13 1 1Again, our ranova analysis shows that the random intercepts by ID_fish are not significant to this model. However, this will be retained because repeated measures were made on the same individuals, and so that the model can be included in the below covariance matrix.
Plot SMR vs. log mass mean predicted level
pred_SMR <- ggpredict(SMR_mod, terms = c("log_mass", "Temp_class"))
ggplot() +
geom_errorbar(data = pred_SMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) + # add error bars showing 95% confidence intervals
geom_point(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group)) +
geom_line(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group)) +
labs(x = "log_mass", y = "Predicted SMR", color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10() 
F test on interaction effect
anova(SMR_mod, type = "3")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log_mass 11.1822 11.1822 1 274.32 1321.4079 < 2.2e-16
## as.factor(Temp_class) 0.6359 0.6359 1 235.74 75.1494 7.134e-16
## Sex 0.0318 0.0159 2 125.25 1.8793 0.156984
## Density 0.0039 0.0039 1 173.19 0.4618 0.497689
## log_mass:as.factor(Temp_class) 0.0861 0.0861 1 277.92 10.1700 0.001591
##
## log_mass ***
## as.factor(Temp_class) ***
## Sex
## Density
## log_mass:as.factor(Temp_class) **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Generate fitted data
z_SMR<-augment(SMR_mod) # get variables and fitted values from model
head(z_SMR)## # A tibble: 6 × 18
## .rownames `log10(SMR)` log_mass `as.factor(Temp_class)` Sex Density ID_fish
## <chr> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
## 1 1 -0.0733 0.423 23 F 8 81
## 2 2 -0.546 -0.0915 23 M 8 82
## 3 3 -0.276 0.250 23 F 8 83
## 4 4 -0.0948 0.436 23 M 8 84
## 5 5 -0.151 0.360 23 M 8 85
## 6 7 -0.296 0.310 23 M 8 87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## # .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## # .weights <dbl>, .wtres <dbl>z_SMR <- z_SMR %>%
mutate(
pwr.fitted = 10^(.fitted)
)
##### check assumptions: plot all residual values against the fitted
ggplot(z_SMR, aes(x = .fitted, y = .resid)) + # same plot as before
geom_point(size = 2) 
ggplot(z_SMR, aes(x = log_mass, y = SMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`), # note this will look odd if there are additional fixed effects
alpha = 0.6) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10()+
labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))
At this stage, we have found a strong model for SMR which captures variance by treatment and uses centered mass. The next step is to repeat this modeling approach for RMR and MMR.
RMR MODEL
# model using log10 RMR and the interaction of log10 mass and temperature treatment
RMR_mod <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 |ID_fish),
data=ds1, REML=T)
plot(RMR_mod)
summary(RMR_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -1011.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5876 -0.6152 0.0310 0.6508 2.7643
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.001429 0.03781
## Residual 0.006328 0.07955
## Number of obs: 509, groups: ID_fish, 113
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.446375 0.032520 245.248665 -13.726
## log_mass 0.714771 0.025398 410.246000 28.142
## as.factor(Temp_class)23 0.152471 0.021214 293.807345 7.187
## SexF -0.017279 0.015909 157.912118 -1.086
## SexM -0.016549 0.016440 144.900150 -1.007
## Density -0.003568 0.003490 252.616585 -1.022
## log_mass:as.factor(Temp_class)23 -0.138010 0.039679 393.028274 -3.478
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 5.48e-12 ***
## SexF 0.279089
## SexM 0.315800
## Density 0.307526
## log_mass:as.factor(Temp_class)23 0.000561 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.449
## as.fc(T_)23 -0.360 0.595
## SexF -0.385 -0.183 -0.073
## SexM -0.396 -0.067 -0.036 0.756
## Density -0.843 0.155 0.123 0.106 0.095
## lg_:.(T_)23 0.321 -0.630 -0.866 0.043 -0.033 -0.108ranova(RMR_mod)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 9 505.67 -993.33
## (1 | ID_fish) 8 494.08 -972.16 23.174 1 1.48e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
RMR_mod_month <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),
data=ds1, REML=T)
plot(RMR_mod_month) # generate model output
summary(RMR_mod_month) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 + Month_int | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -1011.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6009 -0.6090 0.0215 0.6603 2.7544
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 1.876e-03 0.043312
## Month_int 3.614e-06 0.001901 -1.00
## Residual 6.317e-03 0.079482
## Number of obs: 509, groups: ID_fish, 113
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.445366 0.032513 245.812680 -13.698
## log_mass 0.716091 0.025404 410.454083 28.188
## as.factor(Temp_class)23 0.150874 0.021468 219.587640 7.028
## SexF -0.017492 0.016095 137.340763 -1.087
## SexM -0.017322 0.016597 133.136674 -1.044
## Density -0.003705 0.003468 238.248975 -1.068
## log_mass:as.factor(Temp_class)23 -0.137589 0.039653 393.788850 -3.470
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 2.61e-11 ***
## SexF 0.279008
## SexM 0.298544
## Density 0.286492
## log_mass:as.factor(Temp_class)23 0.000578 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.452
## as.fc(T_)23 -0.360 0.596
## SexF -0.392 -0.187 -0.073
## SexM -0.405 -0.070 -0.032 0.763
## Density -0.838 0.156 0.120 0.110 0.099
## lg_:.(T_)23 0.325 -0.630 -0.870 0.045 -0.033 -0.110
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')ranova(RMR_mod_month) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df
## <none> 11 505.81 -989.62
## Month_int in (1 + Month_int | ID_fish) 9 505.67 -993.33 0.28935 2
## Pr(>Chisq)
## <none>
## Month_int in (1 + Month_int | ID_fish) 0.8653# model with random intercepts by temp treatment
RMR_mod_trtint <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),
data=ds1, REML=T)
plot(RMR_mod_trtint)
summary(RMR_mod_trtint)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 | Temp_class)
## Data: ds1
##
## REML criterion at convergence: -988.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9116 -0.6218 0.0056 0.6363 3.1738
##
## Random effects:
## Groups Name Variance Std.Dev.
## Temp_class (Intercept) 0.000744 0.02728
## Residual 0.007684 0.08766
## Number of obs: 509, groups: Temp_class, 2
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -4.465e-01 3.919e-02 1.691e-08 -11.394
## log_mass 7.095e-01 2.365e-02 5.020e+02 30.001
## as.factor(Temp_class)23 1.387e-01 4.296e-02 6.104e-09 3.229
## SexF -1.615e-02 1.336e-02 5.020e+02 -1.209
## SexM -1.567e-02 1.362e-02 5.020e+02 -1.150
## Density -3.326e-03 2.999e-03 5.020e+02 -1.109
## log_mass:as.factor(Temp_class)23 -1.078e-01 3.651e-02 5.020e+02 -2.954
## Pr(>|t|)
## (Intercept) 1.00000
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 1.00000
## SexF 0.22718
## SexM 0.25062
## Density 0.26781
## log_mass:as.factor(Temp_class)23 0.00328 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.323
## as.fc(T_)23 -0.563 0.276
## SexF -0.284 -0.194 -0.029
## SexM -0.300 -0.065 -0.004 0.797
## Density -0.596 0.114 0.052 0.126 0.111
## lg_:.(T_)23 0.247 -0.639 -0.397 0.045 -0.039 -0.105ranova(RMR_mod_trtint)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 9 494.08 -970.16
## (1 | Temp_class) 8 494.08 -972.16 1.1369e-13 1 1# model with random slopes by temp treatment
RMR_mod_trt <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0 + Temp_class|ID_fish),
data=ds1, REML=T)
plot(RMR_mod_trt)
summary(RMR_mod_trt)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (0 + Temp_class | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -1011.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5461 -0.6075 0.0311 0.6356 2.8170
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish Temp_class18 0.001623 0.04028
## Temp_class23 0.001110 0.03331 0.06
## Residual 0.006328 0.07955
## Number of obs: 509, groups: ID_fish, 113
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.445596 0.032521 243.431641 -13.702
## log_mass 0.715638 0.025664 339.122224 27.884
## as.factor(Temp_class)23 0.150977 0.020968 254.500115 7.200
## SexF -0.017989 0.015932 154.494731 -1.129
## SexM -0.017352 0.016446 142.933042 -1.055
## Density -0.003650 0.003464 233.113869 -1.054
## log_mass:as.factor(Temp_class)23 -0.134603 0.039235 317.874051 -3.431
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 6.76e-12 ***
## SexF 0.260576
## SexM 0.293156
## Density 0.293117
## log_mass:as.factor(Temp_class)23 0.000682 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.453
## as.fc(T_)23 -0.375 0.609
## SexF -0.390 -0.175 -0.072
## SexM -0.398 -0.067 -0.037 0.762
## Density -0.839 0.154 0.128 0.111 0.098
## lg_:.(T_)23 0.334 -0.644 -0.869 0.036 -0.035 -0.114
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluesranova(RMR_mod_trt)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df
## <none> 11 505.87 -989.74
## Temp_class in (0 + Temp_class | ID_fish) 9 505.67 -993.33 0.41398 2
## Pr(>Chisq)
## <none>
## Temp_class in (0 + Temp_class | ID_fish) 0.813# plot RMR vs. log mass mean predicted level
pred_RMR <- ggpredict(RMR_mod, terms = c("log_mass", "Temp_class"))
ggplot() +
geom_errorbar(data = pred_RMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) + # add error bars showing 95% confidence intervals
geom_point(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group)) +
geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group)) +
labs(x = "log_mass_centre", y = "Predicted RMR", color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10() 
# plot the raw data as a cloud around the predictions
ggplot() +
geom_point(data = ds1, aes(x = log_mass, y = RMR, group = factor(Temp_class), color = factor(Temp_class)),
alpha = 0.3, position = position_jitter(width = 0.1)) +
geom_point(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group),
size = 3) +
geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group),
size = 2) +
labs(x = "log_mass_centre", y = "Predicted RMR", color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10() 
# F test on interaction effect
anova(RMR_mod, type = "3")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log_mass 6.5282 6.5282 1 386.99 1031.6186 < 2.2e-16
## as.factor(Temp_class) 0.3269 0.3269 1 293.81 51.6568 5.482e-12
## Sex 0.0080 0.0040 2 120.77 0.6298 0.5344197
## Density 0.0066 0.0066 1 252.62 1.0455 0.3075257
## log_mass:as.factor(Temp_class) 0.0766 0.0766 1 393.03 12.0973 0.0005614
##
## log_mass ***
## as.factor(Temp_class) ***
## Sex
## Density
## log_mass:as.factor(Temp_class) ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# generate fitted data
z_RMR<-augment(RMR_mod) # get variables and fitted values from model
head(z_RMR)## # A tibble: 6 × 18
## .rownames `log10(RMR)` log_mass `as.factor(Temp_class)` Sex Density ID_fish
## <chr> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
## 1 1 -0.0212 0.423 23 F 8 81
## 2 2 -0.386 -0.0915 23 M 8 82
## 3 3 -0.147 0.250 23 F 8 83
## 4 4 0.00674 0.436 23 M 8 84
## 5 5 -0.0573 0.360 23 M 8 85
## 6 7 -0.161 0.310 23 M 8 87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## # .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## # .weights <dbl>, .wtres <dbl>z_RMR <- z_RMR %>%
mutate(
pwr.fitted = 10^(.fitted)
)
# check assumptions: plot all residual values against the fitted
ggplot(z_RMR, aes(x = .fitted, y = .resid)) + # same plot as before
geom_point(size = 2) 
ggplot(z_RMR, aes(x = log_mass, y = RMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),
alpha = 0.6) +
theme_bw() +
scale_y_log10() +
scale_color_manual(values = c("#78B7C5", "#F21A00"))
MMR MODEL
# model using log10 MMR and the interaction of log10 mass and temperature treatment
MMR_mod <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 |ID_fish),
data=ds1, REML=T)
plot(MMR_mod)
summary(MMR_mod) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -1166.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5093 -0.5974 0.0398 0.5682 3.4530
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_fish (Intercept) 0.001423 0.03772
## Residual 0.004466 0.06683
## Number of obs: 509, groups: ID_fish, 113
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.013632 0.028910 251.274131 -0.472
## log_mass 0.696766 0.022112 432.420918 31.511
## as.factor(Temp_class)23 0.047231 0.018743 292.436526 2.520
## SexF -0.015674 0.014420 139.112250 -1.087
## SexM 0.008804 0.014953 127.598561 0.589
## Density -0.011717 0.003096 271.720060 -3.784
## log_mass:as.factor(Temp_class)23 -0.084810 0.034606 419.365154 -2.451
## Pr(>|t|)
## (Intercept) 0.63766
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 0.01227 *
## SexF 0.27891
## SexM 0.55705
## Density 0.00019 ***
## log_mass:as.factor(Temp_class)23 0.01466 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.448
## as.fc(T_)23 -0.354 0.585
## SexF -0.387 -0.178 -0.074
## SexM -0.395 -0.066 -0.042 0.747
## Density -0.843 0.166 0.123 0.101 0.090
## lg_:.(T_)23 0.314 -0.628 -0.852 0.041 -0.032 -0.108ranova(MMR_mod)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 9 583.23 -1148.5
## (1 | ID_fish) 8 567.20 -1118.4 32.075 1 1.483e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
MMR_mod_month <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),
data=ds1, REML=T)
plot(MMR_mod_month) # generate model output
summary(MMR_mod_month) # plot model residuals## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 + Month_int | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -1170.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1106 -0.5810 0.0524 0.5567 3.4994
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish (Intercept) 0.0032626 0.05712
## Month_int 0.0001838 0.01356 -0.76
## Residual 0.0040761 0.06384
## Number of obs: 509, groups: ID_fish, 113
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.013991 0.029119 236.454925 -0.480
## log_mass 0.694684 0.022919 329.612730 30.310
## as.factor(Temp_class)23 0.043012 0.019444 201.617723 2.212
## SexF -0.017277 0.014500 130.930308 -1.192
## SexM 0.006459 0.014973 125.441099 0.431
## Density -0.011152 0.003128 218.545598 -3.566
## log_mass:as.factor(Temp_class)23 -0.077658 0.036235 275.968002 -2.143
## Pr(>|t|)
## (Intercept) 0.631343
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 0.028082 *
## SexF 0.235604
## SexM 0.666953
## Density 0.000446 ***
## log_mass:as.factor(Temp_class)23 0.032973 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.433
## as.fc(T_)23 -0.350 0.590
## SexF -0.399 -0.182 -0.067
## SexM -0.409 -0.065 -0.028 0.760
## Density -0.835 0.126 0.104 0.114 0.100
## lg_:.(T_)23 0.311 -0.622 -0.870 0.038 -0.038 -0.091ranova(MMR_mod_month) # anova-like test for significance of model components## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 585.47 -1148.9
## Month_int in (1 + Month_int | ID_fish) 9 583.23 -1148.5 4.47 2 0.107# model with random intercepts by temp treatment
MMR_mod_trtint <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),
data=ds1, REML=T)
plot(MMR_mod_trtint)
summary(MMR_mod_trtint) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (1 | Temp_class)
## Data: ds1
##
## REML criterion at convergence: -1134.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2055 -0.6028 0.0431 0.6245 2.9068
##
## Random effects:
## Groups Name Variance Std.Dev.
## Temp_class (Intercept) 0.000000 0.00000
## Residual 0.005742 0.07578
## Number of obs: 509, groups: Temp_class, 2
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.030099 0.024328 502.000000 -1.237
## log_mass 0.699658 0.020445 502.000000 34.222
## as.factor(Temp_class)23 0.038596 0.016348 502.000000 2.361
## SexF -0.015362 0.011546 502.000000 -1.330
## SexM 0.009194 0.011776 502.000000 0.781
## Density -0.009563 0.002592 502.000000 -3.689
## log_mass:as.factor(Temp_class)23 -0.065391 0.031558 502.000000 -2.072
## Pr(>|t|)
## (Intercept) 0.216586
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 0.018612 *
## SexF 0.183961
## SexM 0.435357
## Density 0.000249 ***
## log_mass:as.factor(Temp_class)23 0.038765 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.449
## as.fc(T_)23 -0.382 0.626
## SexF -0.396 -0.194 -0.065
## SexM -0.418 -0.065 -0.009 0.797
## Density -0.831 0.114 0.119 0.126 0.111
## lg_:.(T_)23 0.344 -0.639 -0.903 0.045 -0.039 -0.105
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')ranova(MMR_mod_trtint)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df Pr(>Chisq)
## <none> 9 567.2 -1116.4
## (1 | Temp_class) 8 567.2 -1118.4 -2.2737e-13 1 1# model with random slopes by temp treatment
MMR_mod_trt <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0 + Temp_class|ID_fish),
data=ds1, REML=T)
plot(MMR_mod_trt)
summary(MMR_mod_trt) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +
## (0 + Temp_class | ID_fish)
## Data: ds1
##
## REML criterion at convergence: -1168.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3634 -0.6080 0.0343 0.5745 3.4158
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID_fish Temp_class18 0.0018262 0.04273
## Temp_class23 0.0008205 0.02865 -0.12
## Residual 0.0044550 0.06675
## Number of obs: 509, groups: ID_fish, 113
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -0.016069 0.028809 246.667733 -0.558
## log_mass 0.698513 0.022560 378.875938 30.962
## as.factor(Temp_class)23 0.045074 0.018239 264.045593 2.471
## SexF -0.017409 0.014412 136.457536 -1.208
## SexM 0.006586 0.014893 126.981343 0.442
## Density -0.011265 0.003029 231.897216 -3.719
## log_mass:as.factor(Temp_class)23 -0.079325 0.033707 333.250899 -2.353
## Pr(>|t|)
## (Intercept) 0.577491
## log_mass < 2e-16 ***
## as.factor(Temp_class)23 0.014095 *
## SexF 0.229152
## SexM 0.659105
## Density 0.000251 ***
## log_mass:as.factor(Temp_class)23 0.019184 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_mss a.(T_) SexF SexM Densty
## log_mass -0.457
## as.fc(T_)23 -0.390 0.615
## SexF -0.400 -0.158 -0.068
## SexM -0.401 -0.065 -0.044 0.762
## Density -0.833 0.162 0.134 0.111 0.098
## lg_:.(T_)23 0.346 -0.660 -0.859 0.024 -0.039 -0.121
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluesranova(MMR_mod_trt)## ANOVA-like table for random-effects: Single term deletions
##
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
## npar logLik AIC LRT Df
## <none> 11 584.32 -1146.7
## Temp_class in (0 + Temp_class | ID_fish) 9 583.23 -1148.5 2.1768 2
## Pr(>Chisq)
## <none>
## Temp_class in (0 + Temp_class | ID_fish) 0.3368# plot MMR vs. log mass mean predicted level
pred_MMR <- ggpredict(MMR_mod, terms = c("log_mass", "Temp_class"))
ggplot() +
geom_errorbar(data = pred_MMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) + # add error bars showing 95% confidence intervals
geom_point(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group)) +
geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group)) +
labs(x = "log_mass", y = "Predicted MMR", color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10() 
# Plot the raw data as a cloud around the predictions
ggplot() +
geom_point(data = ds1, aes(x = log_mass, y = MMR, group = factor(Temp_class), color = factor(Temp_class)),
alpha = 0.3, position = position_jitter(width = 0.1)) +
geom_point(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group),
size = 3) +
geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group),
size = 2) +
labs(x = "log_mass_centre", y = "Predicted MMR", color = "Temp_class") +
theme_bw() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10() 
# F test on interaction effect
anova(MMR_mod, type = "3")## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log_mass 6.2104 6.2104 1 412.94 1390.6566 < 2.2e-16
## as.factor(Temp_class) 0.0284 0.0284 1 292.44 6.3497 0.0122725
## Sex 0.0251 0.0126 2 108.54 2.8111 0.0645219
## Density 0.0640 0.0640 1 271.72 14.3207 0.0001897
## log_mass:as.factor(Temp_class) 0.0268 0.0268 1 419.37 6.0062 0.0146630
##
## log_mass ***
## as.factor(Temp_class) *
## Sex .
## Density ***
## log_mass:as.factor(Temp_class) *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Generate fitted data
z_MMR<-augment(MMR_mod) # get variables and fitted values from model
head(z_MMR)## # A tibble: 6 × 18
## .rownames `log10(MMR)` log_mass `as.factor(Temp_class)` Sex Density ID_fish
## <chr> <dbl> <dbl> <fct> <fct> <dbl> <dbl>
## 1 1 0.154 0.423 23 F 8 81
## 2 2 0.0458 -0.0915 23 M 8 82
## 3 3 0.159 0.250 23 F 8 83
## 4 4 0.215 0.436 23 M 8 84
## 5 5 0.227 0.360 23 M 8 85
## 6 7 0.0919 0.310 23 M 8 87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## # .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## # .weights <dbl>, .wtres <dbl>z_MMR <- z_MMR %>%
mutate(
pwr.fitted = 10^(.fitted)
)
# check assumptions: plot all residual values against the fitted
ggplot(z_MMR, aes(x = .fitted, y = .resid)) + # same plot as before
geom_point(size = 2) 
ggplot(z_MMR, aes(x = log_mass, y = MMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),
alpha = 0.6) +
theme_bw() +
scale_y_log10()+
scale_color_manual(values = c("#78B7C5", "#F21A00"))
FIGURE 2
Generate a figure that displays mass-scaling relationships for all three metabolic rates at both temperatures
### FIG. 2 ###
fig2.1 <- ggplot() +
geom_ribbon(data = pred_SMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) + # add error bars showing 95% confidence intervals
geom_line(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group),
size = 1) +
geom_point(data = ds1, aes(x = log_mass, y = SMR, group = factor(Temp_class), color = factor(Temp_class)),
alpha = 0.25, position = position_jitter(width = 0.1)) +
labs(x = "log10(mass (g))", y = expression(SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
xlim(-0.25, 1.1) +
scale_y_log10(limits = c(0.15, 1.6))
fig2.2 <- ggplot() +
geom_ribbon(data = pred_RMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) + # add error bars showing 95% confidence intervals
geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group),
size = 1) +
geom_point(data = ds1, aes(x = log_mass, y = RMR, group = factor(Temp_class), color = factor(Temp_class)),
alpha = 0.25, position = position_jitter(width = 0.1)) +
labs(x = "log10(mass (g))", y = expression(RMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
xlim(-0.25, 1.1) +
scale_y_log10(limits=c(0.25,2.1))
fig2.3 <- ggplot() +
geom_ribbon(data = pred_MMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) + # add error bars showing 95% confidence intervals
geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group),
size = 1) +
geom_point(data = ds1, aes(x = log_mass, y = MMR, group = factor(Temp_class), color = factor(Temp_class)),
alpha = 0.25, position = position_jitter(width = 0.1)) +
labs(x = expression(log[10]~(mass~(g))), y = expression(MMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
xlim(-0.25, 1.1) +
scale_y_log10(limits=c(0.5, 5))
fig2.4 <- ggplot(z_SMR, aes(x = log_mass, y = SMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`), # note this will look odd if there are additional fixed effects
alpha = 0.6) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10(limits = c(0.15, 1.6))+
xlim(-0.25, 1.1) +
labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))
fig2.5 <- ggplot(z_RMR, aes(x = log_mass, y = RMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`), # note this will look odd if there are additional fixed effects
alpha = 0.6) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10(limits=c(0.25,2.1))+
xlim(-0.25, 1.1) +
labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~RMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))
fig2.6 <- ggplot(z_MMR, aes(x = log_mass, y = MMR, group = ID_fish)) +
geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`), # note this will look odd if there are additional fixed effects
alpha = 0.6) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) +
scale_y_log10(limits=c(0.5, 5))+
xlim(-0.25, 1.1) +
labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~MMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))
library(ggpubr)
fig2 <- ggarrange(print(fig2.1 + rremove("xlab")),
print(fig2.2 + rremove("xlab")),
print(fig2.3 + rremove("xlab")),
print(fig2.4 + rremove("xlab") + rremove("ylab")),
print(fig2.5 + rremove("xlab") + rremove("ylab")),
print(fig2.6 + rremove("xlab") + rremove("ylab")),
ncol = 3, nrow = 2, widths = c(5, 5, 5),
align = "v",
common.legend = TRUE,
legend = "right")
require(grid)
annotate_figure(fig2, bottom = textGrob(expression(log[10]~(mass~(g))), gp = gpar(cex = 1.1)))
setwd(figures_wd)
jpeg(file = "Fig2.jpeg", width = 1440, height=960, units = "px")
annotate_figure(fig2, bottom = textGrob(expression(log[10]~(mass~(g))), gp = gpar(cex = 1.1)))
dev.off()## png
## 2UNIVARIATE MODEL OF FAT MASS
#subset ds1 to only include fish for which fat mass was measured
ds_fat <- ds1[which(!is.na(ds1$log_fatmass)),]
fat_mod <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + Sex + Density, data=ds_fat)
fat_mod1 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + SMR:Temp_class + Sex + Density, data=ds_fat)
fat_mod2 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + RMR:Temp_class + Sex + Density, data=ds_fat)
fat_mod3 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + MMR:Temp_class + Sex + Density, data=ds_fat)
summary(fat_mod) # generate model output##
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class +
## Sex + Density, data = ds_fat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.83753 -0.09771 0.03729 0.12629 0.31222
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.714798 0.064714 -26.498 < 2e-16 ***
## Temp_class23 -0.337119 0.062082 -5.430 9.31e-08 ***
## SexM -0.013164 0.018726 -0.703 0.482
## Density -0.006225 0.006450 -0.965 0.335
## Temp_class18:log_finalmass_centre 1.627573 0.065734 24.760 < 2e-16 ***
## Temp_class23:log_finalmass_centre 2.421200 0.083660 28.941 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1826 on 442 degrees of freedom
## (115 observations deleted due to missingness)
## Multiple R-squared: 0.7806, Adjusted R-squared: 0.7781
## F-statistic: 314.4 on 5 and 442 DF, p-value: < 2.2e-16plot(fat_mod) # residuals plot
summary(fat_mod1) # generate model output##
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class +
## SMR:Temp_class + Sex + Density, data = ds_fat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.82513 -0.08758 0.03480 0.12476 0.34453
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.709110 0.064990 -26.298 < 2e-16 ***
## Temp_class23 -0.324979 0.062354 -5.212 2.88e-07 ***
## SexM -0.015560 0.018727 -0.831 0.4065
## Density -0.006747 0.006456 -1.045 0.2966
## Temp_class18:log_finalmass_centre 1.617051 0.075888 21.308 < 2e-16 ***
## Temp_class23:log_finalmass_centre 2.574283 0.114323 22.518 < 2e-16 ***
## Temp_class18:SMR 0.008286 0.043880 0.189 0.8503
## Temp_class23:SMR -0.136993 0.069860 -1.961 0.0505 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1822 on 440 degrees of freedom
## (115 observations deleted due to missingness)
## Multiple R-squared: 0.7825, Adjusted R-squared: 0.779
## F-statistic: 226.1 on 7 and 440 DF, p-value: < 2.2e-16plot(fat_mod) # residuals plot
summary(fat_mod2) # generate model output##
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class +
## RMR:Temp_class + Sex + Density, data = ds_fat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.82329 -0.08922 0.03455 0.12615 0.32535
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.7069333 0.0655146 -26.054 < 2e-16 ***
## Temp_class23 -0.3273298 0.0627418 -5.217 2.8e-07 ***
## SexM -0.0162856 0.0188321 -0.865 0.388
## Density -0.0068568 0.0064714 -1.060 0.290
## Temp_class18:log_finalmass_centre 1.6226214 0.0766152 21.179 < 2e-16 ***
## Temp_class23:log_finalmass_centre 2.5482306 0.1180100 21.593 < 2e-16 ***
## Temp_class18:RMR 0.0006086 0.0348560 0.017 0.986
## Temp_class23:RMR -0.0865418 0.0567139 -1.526 0.128
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1825 on 440 degrees of freedom
## (115 observations deleted due to missingness)
## Multiple R-squared: 0.7817, Adjusted R-squared: 0.7782
## F-statistic: 225.1 on 7 and 440 DF, p-value: < 2.2e-16plot(fat_mod) # residuals plot
summary(fat_mod3) # generate model output##
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class +
## MMR:Temp_class + Sex + Density, data = ds_fat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.83208 -0.08588 0.03875 0.11949 0.32539
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.694523 0.066308 -25.555 < 2e-16 ***
## Temp_class23 -0.338613 0.063490 -5.333 1.55e-07 ***
## SexM -0.013048 0.018730 -0.697 0.486
## Density -0.007680 0.006523 -1.177 0.240
## Temp_class18:log_finalmass_centre 1.679581 0.079663 21.084 < 2e-16 ***
## Temp_class23:log_finalmass_centre 2.489030 0.110451 22.535 < 2e-16 ***
## Temp_class18:MMR -0.020546 0.017446 -1.178 0.240
## Temp_class23:MMR -0.024615 0.026652 -0.924 0.356
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1825 on 440 degrees of freedom
## (115 observations deleted due to missingness)
## Multiple R-squared: 0.7817, Adjusted R-squared: 0.7782
## F-statistic: 225 on 7 and 440 DF, p-value: < 2.2e-16plot(fat_mod) # residuals plot
pred_fat <- ggpredict(fat_mod, terms = c("log_finalmass_centre", "Temp_class"))
# Plot the raw data
ggplot() +
geom_point(data = ds_fat, aes(x = log_finalmass_centre, y = log_fatmass, group = factor(Temp_class), color = factor(Temp_class))) +
geom_line(data = pred_fat, aes(x = x, y = predicted, color = group, group = group)) +
labs(x = expression(log[10]~(centred~mass~(g))), y = expression(log[10]~(fat~mass~(g))), color = expression("Temperature ("*degree*C*")")) +
theme_classic() +
scale_color_manual(values = c("#78B7C5", "#F21A00")) 
# F test on interaction term
anova(fat_mod)## Analysis of Variance Table
##
## Response: log_fatmass
## Df Sum Sq Mean Sq F value Pr(>F)
## Temp_class 1 0.214 0.2144 6.4335 0.01154 *
## Sex 1 1.014 1.0144 30.4389 5.877e-08 ***
## Density 1 0.009 0.0086 0.2594 0.61080
## Temp_class:log_finalmass_centre 2 51.159 25.5793 767.5449 < 2.2e-16 ***
## Residuals 442 14.730 0.0333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fat_mod1)## Analysis of Variance Table
##
## Response: log_fatmass
## Df Sum Sq Mean Sq F value Pr(>F)
## Temp_class 1 0.214 0.2144 6.4609 0.01137 *
## Sex 1 1.014 1.0144 30.5685 5.534e-08 ***
## Density 1 0.009 0.0086 0.2605 0.61004
## Temp_class:log_finalmass_centre 2 51.159 25.5793 770.8136 < 2.2e-16 ***
## Temp_class:SMR 2 0.129 0.0644 1.9412 0.14476
## Residuals 440 14.601 0.0332
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fat_mod2)## Analysis of Variance Table
##
## Response: log_fatmass
## Df Sum Sq Mean Sq F value Pr(>F)
## Temp_class 1 0.214 0.2144 6.4383 0.01151 *
## Sex 1 1.014 1.0144 30.4615 5.827e-08 ***
## Density 1 0.009 0.0086 0.2596 0.61067
## Temp_class:log_finalmass_centre 2 51.159 25.5793 768.1164 < 2.2e-16 ***
## Temp_class:RMR 2 0.078 0.0388 1.1646 0.31302
## Residuals 440 14.653 0.0333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fat_mod3)## Analysis of Variance Table
##
## Response: log_fatmass
## Df Sum Sq Mean Sq F value Pr(>F)
## Temp_class 1 0.214 0.2144 6.4367 0.01152 *
## Sex 1 1.014 1.0144 30.4540 5.848e-08 ***
## Density 1 0.009 0.0086 0.2595 0.61071
## Temp_class:log_finalmass_centre 2 51.159 25.5793 767.9250 < 2.2e-16 ***
## Temp_class:MMR 2 0.074 0.0370 1.1094 0.33066
## Residuals 440 14.656 0.0333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This model is excluded from the following covariance matrix because we were only able to collect one measurement of fat mass from each fish, all at the very end of the experiment. Without an initial measurement of fat mass (logistically impossible as it was a lethal measurement), we are unable to assess within-subject residual correlation. For this reason, we have assessed it separately.
COVARIANCE MATRIX
build covariance matrix with brms
RN <- bf(mass1 ~ 0 + Temp_class + Month_centred_factor:Temp_class + Sex + Density + (1 + Month_centred |a| ID_fish)) +
bf(SMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) +
bf(RMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) +
bf(MMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) +
set_rescor(T)
prior1 <- c(set_prior("normal(0,2)", class = "b", resp = c("mass1", "SMR1", "RMR1","MMR1")),
set_prior("cauchy(0,2)", class = "sd", resp = c("mass1", "SMR1", "RMR1","MMR1")),
set_prior("lkj(2)", class = "cor"))
fit1 <- brm(formula = RN,
data = ds1,
prior = prior1,
seed = 42, warmup = 500, iter = 10000, chains = 4, cores = 4,
control=list(adapt_delta=0.97, max_treedepth = 15),
save_ranef = T
)
print(fit1, digits = 3)## Family: MV(gaussian, gaussian, gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: mass1 ~ 0 + Temp_class + Month_centred_factor:Temp_class + Sex + Density + (1 + Month_centred | a | ID_fish)
## SMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish)
## RMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish)
## MMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish)
## Data: ds1 (Number of observations: 509)
## Draws: 4 chains, each with iter = 10000; warmup = 500; thin = 1;
## total post-warmup draws = 38000
##
## Multilevel Hyperparameters:
## ~ID_fish (Number of levels: 113)
## Estimate Est.Error l-95% CI u-95% CI
## sd(mass1_Intercept) 0.180 0.013 0.156 0.207
## sd(mass1_Month_centred) 0.042 0.004 0.035 0.051
## sd(SMR1_Intercept) 0.015 0.007 0.002 0.028
## sd(RMR1_Intercept) 0.033 0.005 0.024 0.043
## sd(MMR1_Intercept) 0.037 0.005 0.027 0.047
## cor(mass1_Intercept,mass1_Month_centred) -0.350 0.097 -0.529 -0.147
## cor(mass1_Intercept,SMR1_Intercept) -0.300 0.273 -0.758 0.300
## cor(mass1_Month_centred,SMR1_Intercept) 0.434 0.243 -0.145 0.819
## cor(mass1_Intercept,RMR1_Intercept) -0.055 0.154 -0.355 0.245
## cor(mass1_Month_centred,RMR1_Intercept) 0.252 0.142 -0.034 0.518
## cor(SMR1_Intercept,RMR1_Intercept) 0.512 0.256 -0.143 0.852
## cor(mass1_Intercept,MMR1_Intercept) -0.090 0.150 -0.375 0.209
## cor(mass1_Month_centred,MMR1_Intercept) 0.311 0.135 0.034 0.561
## cor(SMR1_Intercept,MMR1_Intercept) 0.169 0.278 -0.416 0.669
## cor(RMR1_Intercept,MMR1_Intercept) 0.416 0.146 0.111 0.681
## Rhat Bulk_ESS Tail_ESS
## sd(mass1_Intercept) 1.000 6949 13169
## sd(mass1_Month_centred) 1.000 7814 16260
## sd(SMR1_Intercept) 1.001 5973 8296
## sd(RMR1_Intercept) 1.002 2822 6513
## sd(MMR1_Intercept) 1.000 13058 18056
## cor(mass1_Intercept,mass1_Month_centred) 1.001 7166 14793
## cor(mass1_Intercept,SMR1_Intercept) 1.000 17582 20725
## cor(mass1_Month_centred,SMR1_Intercept) 1.000 20046 15852
## cor(mass1_Intercept,RMR1_Intercept) 1.000 14263 21929
## cor(mass1_Month_centred,RMR1_Intercept) 1.000 11969 21026
## cor(SMR1_Intercept,RMR1_Intercept) 1.003 1489 2408
## cor(mass1_Intercept,MMR1_Intercept) 1.000 18123 25293
## cor(mass1_Month_centred,MMR1_Intercept) 1.000 14145 21953
## cor(SMR1_Intercept,MMR1_Intercept) 1.001 1933 3914
## cor(RMR1_Intercept,MMR1_Intercept) 1.000 7547 17439
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI
## mass1_Temp_class18 0.175 0.061 0.056 0.294
## mass1_Temp_class23 0.108 0.062 -0.014 0.231
## mass1_SexF 0.121 0.045 0.033 0.208
## mass1_SexM 0.050 0.048 -0.043 0.143
## mass1_Density 0.010 0.006 -0.001 0.022
## mass1_Temp_class18:Month_centred_factor1 0.160 0.012 0.137 0.183
## mass1_Temp_class23:Month_centred_factor1 0.128 0.015 0.098 0.157
## mass1_Temp_class18:Month_centred_factor2 0.224 0.015 0.194 0.255
## mass1_Temp_class23:Month_centred_factor2 0.217 0.019 0.180 0.255
## mass1_Temp_class18:Month_centred_factor3 0.168 0.021 0.126 0.209
## mass1_Temp_class23:Month_centred_factor3 0.225 0.025 0.176 0.273
## mass1_Temp_class18:Month_centred_factor4 0.378 0.025 0.329 0.428
## mass1_Temp_class23:Month_centred_factor4 0.355 0.031 0.295 0.415
## SMR1_Temp_class18 -0.845 0.036 -0.917 -0.775
## SMR1_Temp_class23 -0.635 0.035 -0.704 -0.566
## SMR1_SexF -0.028 0.015 -0.057 0.001
## SMR1_SexM -0.015 0.015 -0.044 0.015
## SMR1_Density -0.002 0.003 -0.008 0.005
## SMR1_Temp_class18:log_mass_centre 0.803 0.031 0.743 0.865
## SMR1_Temp_class23:log_mass_centre 0.675 0.035 0.607 0.746
## RMR1_Temp_class18 -0.607 0.035 -0.677 -0.538
## RMR1_Temp_class23 -0.427 0.034 -0.495 -0.360
## RMR1_SexF -0.017 0.015 -0.047 0.013
## RMR1_SexM -0.016 0.016 -0.047 0.014
## RMR1_Density -0.005 0.003 -0.011 0.002
## RMR1_Temp_class18:log_mass_centre 0.714 0.029 0.657 0.770
## RMR1_Temp_class23:log_mass_centre 0.584 0.033 0.520 0.649
## MMR1_Temp_class18 -0.188 0.033 -0.253 -0.122
## MMR1_Temp_class23 -0.120 0.033 -0.185 -0.056
## MMR1_SexF -0.017 0.014 -0.045 0.011
## MMR1_SexM 0.008 0.015 -0.021 0.037
## MMR1_Density -0.011 0.003 -0.017 -0.005
## MMR1_Temp_class18:log_mass_centre 0.705 0.027 0.652 0.758
## MMR1_Temp_class23:log_mass_centre 0.622 0.031 0.561 0.682
## Rhat Bulk_ESS Tail_ESS
## mass1_Temp_class18 1.000 5589 11325
## mass1_Temp_class23 1.001 5162 9990
## mass1_SexF 1.001 3826 7176
## mass1_SexM 1.001 3950 7862
## mass1_Density 1.000 11273 17979
## mass1_Temp_class18:Month_centred_factor1 1.000 12484 22473
## mass1_Temp_class23:Month_centred_factor1 1.000 13650 23584
## mass1_Temp_class18:Month_centred_factor2 1.001 7445 15542
## mass1_Temp_class23:Month_centred_factor2 1.000 7430 15925
## mass1_Temp_class18:Month_centred_factor3 1.001 6026 12963
## mass1_Temp_class23:Month_centred_factor3 1.001 6060 13153
## mass1_Temp_class18:Month_centred_factor4 1.001 5275 10983
## mass1_Temp_class23:Month_centred_factor4 1.001 5588 11585
## SMR1_Temp_class18 1.000 11492 19639
## SMR1_Temp_class23 1.000 11580 20326
## SMR1_SexF 1.000 15446 23010
## SMR1_SexM 1.000 16118 24116
## SMR1_Density 1.000 18842 26708
## SMR1_Temp_class18:log_mass_centre 1.000 12958 21232
## SMR1_Temp_class23:log_mass_centre 1.000 14744 22700
## RMR1_Temp_class18 1.000 10436 18121
## RMR1_Temp_class23 1.001 10171 17456
## RMR1_SexF 1.000 13066 21330
## RMR1_SexM 1.000 13478 21505
## RMR1_Density 1.000 16399 24852
## RMR1_Temp_class18:log_mass_centre 1.000 12740 20892
## RMR1_Temp_class23:log_mass_centre 1.000 14206 21255
## MMR1_Temp_class18 1.000 14348 22135
## MMR1_Temp_class23 1.000 14599 22742
## MMR1_SexF 1.000 17061 25325
## MMR1_SexM 1.000 17841 24465
## MMR1_Density 1.000 20724 27685
## MMR1_Temp_class18:log_mass_centre 1.000 18451 25115
## MMR1_Temp_class23:log_mass_centre 1.000 21389 26997
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_mass1 0.060 0.003 0.055 0.065 1.000 18474 24645
## sigma_SMR1 0.092 0.003 0.086 0.098 1.000 16337 23784
## sigma_RMR1 0.081 0.003 0.075 0.086 1.000 11070 22175
## sigma_MMR1 0.068 0.002 0.063 0.073 1.000 24222 25735
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(mass1,SMR1) -0.038 0.072 -0.180 0.104 1.000 15109 22918
## rescor(mass1,RMR1) -0.066 0.069 -0.202 0.070 1.000 14908 23218
## rescor(SMR1,RMR1) 0.862 0.013 0.836 0.886 1.000 16327 25735
## rescor(mass1,MMR1) -0.103 0.066 -0.231 0.029 1.000 18390 24837
## rescor(SMR1,MMR1) 0.169 0.047 0.075 0.259 1.000 21566 27134
## rescor(RMR1,MMR1) 0.238 0.047 0.145 0.327 1.000 22986 27361
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(fit1)
conditional_effects(fit1)
REPEATABILITY ESTIMATES
Use posterior distribution to generate repeatability estimates
View(posterior::summarise_draws(as_draws_array(fit1))) # view all predicted values of random effects
# random effects are fitted as standard dev, not variance, so we need to square them first
# var names to insert below: SMR, RMR, MMR, growth
# SMR repeatability
Vint_SMR <- (as_draws_array(fit1, variable = "sd_ID_fish__SMR1_Intercept"))^2
Vresid_SMR <- (as_draws_array(fit1, variable = "sigma_SMR1"))^2
R_SMR <- Vint_SMR / (Vint_SMR + Vresid_SMR)
"SMR"; mean(R_SMR); quantile(R_SMR, probs = c(0.025, 0.975))## [1] "SMR"## [1] 0.03092714## 2.5% 97.5%
## 0.0004320659 0.0890740509# RMR repeatability
Vint_RMR <- (as_draws_array(fit1, variable = "sd_ID_fish__RMR1_Intercept"))^2
Vresid_RMR <- (as_draws_array(fit1, variable = "sigma_RMR1"))^2
R_RMR <- Vint_RMR / (Vint_RMR + Vresid_RMR)
"RMR"; mean(R_RMR); quantile(R_RMR, probs = c(0.025, 0.975))## [1] "RMR"## [1] 0.1458303## 2.5% 97.5%
## 0.07860133 0.23501854# MMR repeatability
Vint_MMR <- (as_draws_array(fit1, variable = "sd_ID_fish__MMR1_Intercept"))^2
Vresid_MMR <- (as_draws_array(fit1, variable = "sigma_MMR1"))^2
R_MMR <- Vint_MMR / (Vint_MMR + Vresid_MMR)
"MMR"; mean(R_MMR); quantile(R_MMR, probs = c(0.025, 0.975))## [1] "MMR"## [1] 0.2306318## 2.5% 97.5%
## 0.1311474 0.3357818# Growth rate repeatability
Vint_growth <- (as_draws_array(fit1, variable = "sd_ID_fish__mass1_Month_centred"))^2
Vslope_growth <- (as_draws_array(fit1, variable = "sd_ID_fish__mass1_Intercept"))^2
corr_growth <- (as_draws_array(fit1, variable = "cor_ID_fish__mass1_Intercept__mass1_Month_centred"))
Vresid_growth <- (as_draws_array(fit1, variable = "sigma_mass1"))^2
COVis_growth <- corr_growth*sqrt(Vslope_growth)*sqrt(Vint_growth)
x<- 0 #(can be 0 - 4)
VAR0 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R0_growth <- VAR0/(VAR0 + Vresid_growth)
"Growth - May"; mean(R0_growth); quantile(R0_growth, probs = c(0.025, 0.975))## [1] "Growth - May"## [1] 0.3336915## 2.5% 97.5%
## 0.2450998 0.4309421x<- 1 #(can be 0 - 4)
VAR1 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R1_growth <- VAR1/(VAR1 + Vresid_growth)
"Growth - June"; mean(R1_growth); quantile(R1_growth, probs = c(0.025, 0.975))## [1] "Growth - June"## [1] 0.8876298## 2.5% 97.5%
## 0.8519065 0.9176981x<- 2 #(can be 0 - 4)
VAR2 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R2_growth <- VAR2/(VAR2 + Vresid_growth)
"Growth - July"; mean(R2_growth); quantile(R2_growth, probs = c(0.025, 0.975))## [1] "Growth - July"## [1] 0.9705673## 2.5% 97.5%
## 0.9599677 0.9790658x<- 3 #(can be 0 - 4)
VAR3 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R3_growth <- VAR3/(VAR3 + Vresid_growth)
"Growth - August"; mean(R3_growth); quantile(R3_growth, probs = c(0.025, 0.975))## [1] "Growth - August"## [1] 0.9869686## 2.5% 97.5%
## 0.9821577 0.9907833x<- 4 #(can be 0 - 4)
VAR4 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R4_growth <- VAR4/(VAR4 + Vresid_growth)
"Growth - September"; mean(R4_growth); quantile(R4_growth, probs = c(0.025, 0.975))## [1] "Growth - September"## [1] 0.9927106## 2.5% 97.5%
## 0.9899955 0.9948520